In [1]:
# Utilities
from __future__ import division
import numpy as np
class MemoizeMutable:
"""
Implement memoization with mutable arguments
"""
def __init__(self, fn):
self.fn = fn
self.memo = {}
def __call__(self, *args, **kwds):
import cPickle
str = cPickle.dumps(args, 1)+cPickle.dumps(kwds, 1)
if not self.memo.has_key(str):
self.memo[str] = self.fn(*args, **kwds)
# else:
# print "returning memoized calculation"
return self.memo[str]
# got from http://stackoverflow.com/questions/1208118/using-numpy-to-build-an-array-of-all-combinations-of-two-arrays
# (faster than other alternatives)
def cartesian(arrays, out=None):
"""
Generate a cartesian product of input arrays.
Parameters
----------
arrays : list of array-like
1-D arrays to form the cartesian product of.
out : ndarray
Array to place the cartesian product in.
Returns
-------
out : ndarray
2-D array of shape (M, len(arrays)) containing cartesian products
formed of input arrays.
Examples
--------
>>> cartesian(([1, 2, 3], [4, 5], [6, 7]))
array([[1, 4, 6],
[1, 4, 7],
[1, 5, 6],
[1, 5, 7],
[2, 4, 6],
[2, 4, 7],
[2, 5, 6],
[2, 5, 7],
[3, 4, 6],
[3, 4, 7],
[3, 5, 6],
[3, 5, 7]])
"""
arrays = [np.asarray(x) for x in arrays]
dtype = arrays[0].dtype
n = np.prod([x.size for x in arrays])
if out is None:
out = np.zeros([n, len(arrays)], dtype=dtype)
m = n / arrays[0].size
out[:,0] = np.repeat(arrays[0], m)
if arrays[1:]:
cartesian(arrays[1:], out=out[0:m,1:])
for j in xrange(1, arrays[0].size):
out[j*m:(j+1)*m,1:] = out[0:m,1:]
return out
In [2]:
from collections import defaultdict, namedtuple
@MemoizeMutable
def fareySequence(N, k=1):
"""
Generate Farey sequence of order N, less than 1/k
"""
# assert type(N) == int, "Order (N) must be an integer"
a, b = 0, 1
c, d = 1, N
seq = [(a,b)]
while c/d <= 1/k:
seq.append((c,d))
tmp = int(math.floor((N+b)/d))
a, b, c, d = c, d, tmp*c-a, tmp*d-b
return seq
@MemoizeMutable
def resonanceSequence(N, k):
"""
Compute resonance sequence
Arguments:
- N (int): Order
- k (int): denominator of the farey frequency resonances are attached to
"""
a, b = 0, 1
c, d = k, N-k
seq = [(a,b)]
while d >= 0:
seq.append((c,d))
tmp = int(math.floor((N+b+a)/(d+c)))
a, b, c, d = c, d, tmp*c-a, tmp*d-b
return seq
def rationalApproximation(points, N, tol=1e-3, lowest_order_only=True):
"""
Arguments:
points: 2D (L x 2) points to approximate
N: max order
"""
L,_ = points.shape
# since this solutions assumes a>0, a 'quick' hack to also obtain solutions
# with a < 0 is to flip the dimensions of the points and explore those
# solutions as well
points = np.vstack((points, np.fliplr(points)))
solutions = defaultdict(set)
sequences = {1: set(fareySequence(1))}
for n in range(2, N+1):
sequences[n] = set(fareySequence(n)) - sequences[n-1]
for h,k in fareySequence(N,1):
if 0 in (h,k):
continue
# print h,k
for x,y in resonanceSequence(N, k):
# avoid 0-solutions
if 0 in (x,y):
continue
norm = np.sqrt(x**2+y**2)
n = np.array([ y/norm, x/norm]) * np.ones_like(points)
n[points[:,0] < h/k, 0] *= -1 # points approaching from the left
# nomenclature inspired in http://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation
ap = np.array([h/k, 0]) - points
apn = np.zeros((1,L))
d = np.zeros_like(points)
apn = np.sum(n*ap, 1, keepdims=True)
d = ap - apn*n
## DON'T RETURN IMMEDIATELY; THERE MIGHT BE OTHER SOLUTIONS OF THE SAME ORDER OR HIGHER
indices, = np.nonzero(np.sqrt(np.sum(d*d,1)) <= tol)
for i in indices:
# print "h/k:", h , "/", k
# print "point:", points[i,:]
if points[i,0] >= h/k:
if i<L:
# print "non-flipped >= h/k"
solutions[i].add((x,-y, h*x/k))
# print i, (x,-y, h*x/k)
elif x*(-y)<0: # only consider solutions where (a,b) have different sign for the "flipped" points (the other solutions should have already been found for the non-flipped points)
# print "flipped >= h/k"
solutions[i-L].add((-y, x, h*x/k))
# print i-L, (-y, x, h*x/k)
else:
if i<L:
# print "non-flipped < h/k"
solutions[i].add((x, y, h*x/k))
# print i, (x, y, h*x/k)
elif x*y>0: # only consider solutions where (a,b) have different sign for the "flipped" points (the other solutions should have already been found for the non-flipped points)
# print "flipped < h/k"
solutions[i-L].add((y, x, h*x/k))
# print i-L, (y, x, h*x/k)
if lowest_order_only:
# removed = 0
for k in solutions:
# keep lowest order solutions only
lowest_order = 2*N
s = set([])
for sol in solutions[k]:
K = abs(sol[0])+abs(sol[1])+abs(sol[2])
if K == lowest_order:
s.add(sol)
elif K < lowest_order:
lowest_order = K
# if len(s) > 0:
# print("point: ({},{}) -> removing {} for {}".format(points[k,0], points[k,1], s, sol))
# removed += len(s)
s = set([sol])
solutions[k] = s
# print("Removed {} solutions".format(removed))
return solutions
@MemoizeMutable
def threeFreqMonomials(fj, fi, allow_self_connect=True, N=5, tol=1e-10, lowest_order_only=False):
"""
Arguments:
fj (np.array_like): frequency vector of the source (j in the paper)
fi (np.array_like): frequency vector of the target (i in the paper)
N: max order
tol: tolerance in the point to line distance calculation
"""
from time import time
st = time()
# use 32bits (64 can take too much memory)
fj = np.array([f for f in fj], dtype=np.float32)
fi = np.array([f for f in fi], dtype=np.float32)
Fj, Fi = len(fj), len(fi)
cart_idx = cartesian((np.arange(Fj),
np.arange(Fj),
np.arange(Fi)))
# we care only when y2 > y1
cart_idx = cart_idx[cart_idx[:,1]>cart_idx[:,0]]
if not allow_self_connect:
cart_idx = cart_idx[(cart_idx[:,0] != cart_idx[:,2]) & (cart_idx[:,1] != cart_idx[:,2])]
# actual frequency triplets
cart = np.vstack((fj[cart_idx[:,0]], fj[cart_idx[:,1]], fi[cart_idx[:,2]])).T
nr, _ = cart_idx.shape
# sort in order to get a*x+b*y=c with 0<x,y<1
sorted_idx = np.argsort(cart, axis=1)
cart.sort()
print("a) Elapsed: {} secs".format(time() - st))
all_points = np.zeros((nr, 2), dtype=np.float32)
all_points[:,0] = cart[:,0] / cart[:,2]
all_points[:,1] = cart[:,1] / cart[:,2]
# del cart
print("b) Elapsed: {} secs".format(time() - st))
redundancy_map = defaultdict(list)
for i in xrange(all_points.shape[0]):
redundancy_map[(all_points[i,0],all_points[i,1])].append(i)
del all_points
print("c) Elapsed: {} secs".format(time() - st))
points = np.array([[a,b] for a,b in redundancy_map])
print("d) Elapsed: {} secs".format(time() - st))
exponents = rationalApproximation(points, N, tol=tol, lowest_order_only=lowest_order_only)
print("e) Elapsed: {} secs".format(time() - st))
monomials = [defaultdict(list) for x in fi]
M = namedtuple('Monomials', ['indices', 'exponents'])
# FIXME: is there a way to reduce the number of nested loops?
for k in exponents:
x, y = points[k,0], points[k,1]
all_points_idx = redundancy_map[(x,y)]
sols = exponents[k]
for a, b, c in sols:
for idx in all_points_idx:
j1, j2, i = (cart_idx[idx, 0], cart_idx[idx, 1], cart_idx[idx, 2])
reordered = (sorted_idx[idx,0], sorted_idx[idx,1], sorted_idx[idx,2])
if reordered == (0,1,2):
n1, n2, d = a, b, c
elif reordered == (0,2,1):
# n1, d, n2 = -a, b, c
n1, n2, d = -a, c, b
elif reordered == (2,0,1):
# d, n1, n2 = a, -b, c
n1, n2, d = -b, c, a
else:
raise Exception("Unimplemented order!")
if d < 0:
n1, n2, d = -n1, -n2, -d
monomials[i]['j1'] += [j1]
monomials[i]['j2'] += [j2+len(fj)] # add offset for fast look-up at run time (will use a flattened array)
monomials[i]['i'] += [i+len(fj)+len(fi)] # add offset for fast look-up at run time (will use a flattened array)
monomials[i]['n1'] += [n1]
monomials[i]['n2'] += [n2]
# monomials[i]['d'] += [d]
monomials[i]['d'] += [d-1] # small optimization to avoid extra subtraction in every step when actually running the network
# print i, (a,b,c), (n1, n2, d)
print("f) Elapsed: {} secs".format(time() - st))
# make them 2d arrays for speed up
for i, m in enumerate(monomials):
I = np.array([m['j1'], m['j2'], m['i']], dtype=np.int).T
E = np.array([m['n1'], m['n2'], m['d']], dtype=np.int).T
monomials[i] = M(indices=I, exponents=E)
print("g) Elapsed: {} secs".format(time() - st))
return monomials
In [5]:
from time import time
fc = 2.0
# num_oscs = 321
num_oscs = 11
f1 = fc * np.logspace(-2.5, 2.5, num_oscs, base=2)
f2 = fc * np.logspace(-2.5, 2.5, num_oscs, base=2)
start = time()
monomials2 = threeFreqMonomials(f1, f2, N=3, tol=1e-10)
print "Run in {:.2f} seconds".format(time()-start)
print monomials
print monomials2
In [ ]: